home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / gnu / fpu881 / src6.zoo / csqrt.c < prev    next >
C/C++ Source or Header  |  1991-09-24  |  3KB  |  115 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    csqrt   complex double precision square root
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    csqrt
  29.  *    machine independent routines
  30.  *    complex functions
  31.  *    math libraries
  32.  *
  33.  *  DESCRIPTION
  34.  *
  35.  *    Computes double precision complex square root of
  36.  *    a double precision complex argument.
  37.  *
  38.  *  USAGE
  39.  *
  40.  *    COMPLEX csqrt (z)
  41.  *    COMPLEX z;
  42.  *
  43.  *  REFERENCES
  44.  *
  45.  *    Fortran 77 user's guide, Digital Equipment Corp. pp B-12
  46.  *
  47.  *  RESTRICTIONS
  48.  *
  49.  *    The relative error in the double precision square root
  50.  *    computation is 10**(-30.1) after three applications
  51.  *    of Heron's iteration for the square root.
  52.  *
  53.  *    However, this assumes exact arithmetic in the iterations
  54.  *    and initial approximation.  Additional errors may occur
  55.  *    due to truncation, rounding, or machine precision limits.
  56.  *    
  57.  *  PROGRAMMER
  58.  *
  59.  *    Fred Fish
  60.  *    Tempe, Az 85281
  61.  *    (602) 966-8871
  62.  *
  63.  *  INTERNALS
  64.  *
  65.  *    Computes complex square root of z = x + j y from:
  66.  *
  67.  *        1.    If z = 0 + j 0 then return z.
  68.  *
  69.  *        2.    root = sqrt((dabs(x) + cabs(z)) / 2)
  70.  *
  71.  *        3.    q = y / (2 * root)
  72.  *
  73.  *        4.    If x >= 0 then
  74.  *            csqrt(z) = (root,q)
  75.  *
  76.  *        5.    If x < 0 and y >= 0 then
  77.  *            csqrt(z) = (q,root)
  78.  *
  79.  *        6.    If x < 0 and y < 0 then
  80.  *            csqrt(z) = (-q,root)
  81.  *
  82.  */
  83.  
  84. #include <stdio.h>
  85. #include <pmluser.h>
  86. #include "pml.h"
  87.  
  88.  
  89. COMPLEX csqrt (z)
  90. COMPLEX z;
  91. {
  92.     double root, q;
  93.     extern double dabs(), sqrt(), cabs ();
  94.  
  95.     ENTER ("csqrt");
  96.     DEBUG4 ("csqrtin", "arg %le %le", z.real, z.imag);
  97.     if (z.real != 0.0 || z.imag != 0.0) {
  98.         root = sqrt (0.5 * (dabs (z.real) + cabs (z)));
  99.         q = z.imag / (2.0 * root);
  100.         if (z.real >= 0.0) {
  101.             z.real = root;
  102.         z.imag = q;
  103.         } else if (z.imag < 0.0) {
  104.         z.real = -q;
  105.         z.imag = -root;
  106.         } else {
  107.         z.real = q;
  108.         z.imag = root;
  109.         }
  110.     }
  111.     DEBUG4 ("csqrtout", "result %le %le", z.real, z.imag);
  112.     LEAVE ();
  113.     return (z);
  114. }
  115.